VARIABLE ELIMINATION IN CHEMICAL REACTION NETWORKS 
WITH MASS ACTION KINETICS 



Abstract. We consider chemical reaction networks taken with mass action kinetics. The 
steady states of such a system are solutions to a system of polynomial equations. Even 
for small systems the task of finding the solutions is daunting. We develop an algebraic 
framework and procedure for linear elimination of variables. The procedure reduces the 
variables in the system to a set of "core" variables by eliminating variables corresponding 
to a set of non-interacting species. The steady states are parameterized algebraically by 
the core variables, and a graphical condition is given for when a steady state with positive 
core variables necessarily have all variables positive. Further, we characterize graphically 
the sets of eliminated variables that are constrained by a conservation law and show that 
this conservation law takes a specific form. 

Keywords: semiflow, species graph, non-interacting species, spanning tree, polynomial 
equations 

1. Introduction 

The goal of this work is to discuss hnear ehmination of variables at steady state in Chem- 
ical Reaction Networks (CRNs) taken with mass action kinetics. We use the formalism of 
Chemical Reaction Network Theory (CRNT) that puts CRNs into a mathematical, in par- 
ticular algebraic, framework. CRNT was developed around 40 years ago, mainly by Horn, 
Jackson and Feinberg [HI [9l [12l [T7] . Its usefulness for analysis of CRNs is continuously 
being supported [21 HI [22]. 

We introduce the elimination procedure by going through a specific example. Consider 
five chemical species j4, C, D, E that interact according to the reactions: 

ri: A + 2B ^ D r2:D^A + C : C + D ^ E r^: E A + B. 

The molar concentration of a species S at time t is denoted by cs = cs{t). By employing 
the common assumption that reaction rates are of mass action type, the concentrations 
change with time according to a system of ordinary differential equations (DDEs): 

CA = -kiCAC% + k2CD + k^CE c'b = -2kiCAC% + ^4C£; c'c = k2CD - k^cccn 

c'd = kiCAC% ~ ^'iCD - k'iCcCD c'e = ksCcCD - kiCE, 

where ki denotes the positive rate constant of reaction rj. Observe that ca + c'd + c'e = 0, 
which implies that ca + ce> + ce is constant over time, and fixed by the sum of initial 
concentrations cq = c^(0) -|- cd{0) + ce{0)- The equation cq = ca + cd + ce is called a 
conservation law. 
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We are interested in the steady state solutions of this system, in particular, the positive 
steady state solutions, that is, the solutions for which all concentrations are positive. The 
steady state solutions are found by setting the ODEs to zero. Consider the equations 
c'a = 0, c'd = 0, and c'e = 0: 

(1.1) = -kiCACB + k2CD + hcE, Q = kiCAC% - k2CD - k^CcCD, 

= k^ccCD - kiCE- 

The equations form a system of polynomial equations in caiCb-,cc-,cd-,ce with real coef- 
ficients. None of the equations contain a monomial with more than one of the variables 
CA, cd,ce- Further, the degree of ca, cd, ce is one in all equations. In other words, if we let 
Con = {ki,k2,k3,k4}, then ( 1.1 ) is a linear system of equations in the variables ca,cd,ce 
with coefficients in the field M(ConU{cB, cc}): 

-kic% k2 ki \ 1 ca\ 

kic\ -k2 - kscc cz) = 0. 

k3Cc -k4: J \ CE / 

The column sums of the 3x3 matrix are zero, because of the conserved amount cq. In 
fact, the matrix has rank 2 in ]R(Con U{cb, cc}) and the solutions of the system form a 
line parameterized for example hy cd- 

CA = -, —rCD, CE = CD- 

k2 + k3 k^cc 
These solutions are well-defined in the field M(ConU{cs, cc})- If positive values of c^, cc, cd 
are given, then the steady state values of ca,ce are positive and completely determined. 
Further, since ca + cd + ce = cq, we find that 

CD = CO 1 + + 



k2 + k3 k^cc, 

and conclude that the positive steady states are fully determined by the positive steady 
state solutions of cb,cc- These solutions are found from the equations c'b = and c'c = 
in cb,cc, by substituting the values of ca,cd,ce- Note that the equations can always be 
rewritten in polynomial form. 

Let us now start from the equations c'c = and c'e = 0: k2CD — k^ccCD = and 
ksCcCD — k^CE = 0. This system is linear in the variables cc,ce with coefficients in the 
field M(ConU{c£)}). Further, the system has maximal rank in M(ConU{cD}) and thus has 
a unique solution cc = k2/k3, ce = k2CD/ki in M(ConU{cD}). As above, the variables 
cc,ce can be eliminated and recovered from any positive steady state solution ca,cb,cd 
of the remaining equations. 

The approaches that are used to eliminate the variables in {ca, cd, ce} and {cc,ce} 
differ: In the former case the system is homogeneous, does not have maximal rank, and the 
conservation law is required for full elimination. In the latter the system has maximal rank 
but is not homogeneous. At this point we might ask: What are the similarities between 
the two sets of variables that enable their elimination from the steady state equations? 
What are the differences that lead to different approaches? The species in both sets do 
not interact with each other, that is, they do not appear on the same side of a reaction, 
and further all concentrations have degree one in all the equations in which they appear. 
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Such sets are called non-interacting. However, in the first case the sum of concentrations is 
conserved, and the set is what we call a cut, while in the second case it is not. Importantly, 
the eliminated variables are non-negative whenever the non-eliminated variables (the core 
variables) are positive. Furthermore, the cases in which zero concentrations of the elimi- 
nated variables can occur can be completely characterized. The concentration cb cannot 
be linearly eliminated because it has degree 2 in some equations. 

This reaction system is small compared to real biochemical systems and can be ma- 
nipulated manually. For an arbitrary CRN, most non-interacting sets can be eliminated 
using one of the approaches outlined above, depending on the presence or absence of a 
conserved amount. After reduction, the (positive) steady states are the solutions to poly- 
nomial equations depending on the core variables only. Thus the steady states form an 
algebraic variety in the core variables. 

In this manuscript we discuss a general procedure for linear elimination of variables, 
embracing the two approaches described above. The design of the procedure relies on (a 
specific version of) the species graph. Subsets of species that can be eliminated are non- 
interacting. These sets correspond to a specific type of subgraphs and in any such set, the 
corresponding subgraph encodes the presence or absence of a conservation law relating the 
concentrations in the set. We study the interplay between non-interacting sets, subgraphs, 
and conservation laws and relate subgraph connectedness to minimality of conservation 
laws and the existence of conservation laws to the so-called full subgraphs. Thus, the 
results obtained here are of interest in their own. The Matrix- Tree theorem [2lj is key to 
study positivity of solutions [23] . 

The elimination procedure has interesting potential applications. First of all, essential 
information about the system at steady state is contained in the equations for the core 
variables. Thus, experimental knowledge about the concentrations of the core variables is 
sufficient to explore the system at steady state. Further, the species graph can be used in 
experimental planning by choosing (if possible) a subgraph that optimizes the information 
in the experiment. 

Secondly, two different reaction systems involving the same chemical species can be 
discriminated based on the core variables alone and thus used for model selection. This 
is possible, irrespectively whether the rate constants are known or not |15l I19j . if the two 
algebraic varieties described by the core variables take different forms. A series of measure- 
ments with diff^erent initial concentrations can determine which variety the measurements 
belong to. 

Finally, another potential application concerns the emergence of multiple positive steady 
states in a specific system. Many mathematical tools for detecting whether a system has 
at most one positive steady state exist [H [2l [10]. However, when these fail, it is not 
straightforward to conclude that the system admits more than one positive steady state and 
rate constants need to be found for which this is true. This is typically done by performing 
a random parameter search. Elimination of variables might reduce the computational 
burden substantially and decrease the likelihood of numerical errors. 

This work builds on our previous work on variable elimination in so-called Post- Transla- 
tional Modification (PTM) systems |13j . PTM systems form a special type of biochemical 
reaction networks that are particularly abundant in cell signaling and have been the focus 
of much theoretical research [161 [El ED]. A subclass of PTM systems was studied by 
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Thomson and Gunawardena [23J. The present work extends the ehmination procedure for 
PTM systems to arbitrary CRNs. By doing so, some particularities of PTM systems are 
uncovered to be irrelevant. 

The outline of the paper is the following. We introduce the notation, some preliminaries, 
and CRNs together with their associated mass action ODEs. We proceed to discuss conser- 
vation laws arising from so-called semiflows, with special attention to minimal semiflows. 
Next, the species graph and its relevant subgraphs (full and non-interacting) are defined, 
and we proceed to discuss relations between the subgraphs and semiflows. We then present 
the variable elimination procedure and the reduction of the steady state equations to a 
polynomial system in the core variables. Using the graphical representation, we show that 
positive solutions of the core variables in the reduced system correspond to non-negative 
steady states of the CRN, in which only the eliminated variables can possibly be zero. 



2. Notation 

Let M+ denote the set of positive real numbers (without zero) and M+ the set of non- 
negative real numbers (with zero). Given a finite set £, let be the real vector space of 
formal sums v = X^sef '^eE, with A^; G M. If A^; G M+ (resp. M+) for all E G £, then we 
write V G (resp. M^). 

S-positivity. Let denote the ring of real polynomials in £. A monomial is a 
polynomial of the form ^YIe^s some A G M \ {0} and n^; G No (the natural 

numbers including zero). A non-zero polynomial in M[£^] with non- negative coefficients is 
called S-positive. Any assignment a: £ ^ M4. induces an evaluation map Ca'- M[iS] — )• M. If 
p G M[i?] is S-positive, then ea{p) > 0. 

A rational function / in iS is S-positive if it is a quotient of two S-positive polynomials 
in £. Then Caif) is well-defined and positive for any assignment a: £ ^ M+. In general, a 
rational function / = p/q in zi, Zs and coefficients in M{£) is S-positive if the coefficients 
of p and q are S-positive rational functions in £. Then Caif) is an S-positive rational 
function in M(zi, . . . , Zg), for any assignment a : £ ^ M4.. Assume that E = g{£) for some 
rational function g m. £ = £ \ {E} and E G £. Then, if / is a rational function in £, 
substituting g into / gives / as a rational function in £. 

Graphs and the Matrix- Tree theorem. Let G be a directed graph with node 
set M. A spanning tree r of G is a directed subgraph with node set M and such that 
the corresponding undirected graph is connected and acyclic. Self-loops are by definition 
excluded from a spanning tree. There is a (unique) undirected path between any two nodes 
in a spanning tree [6]. We say that the spanning tree r is rooted at a node v if the unique 
path between any node w and v is directed from w to v. As a consequence, v is the only 
node in r with no edges of the form v ^ w (called out-edges). Further, there is no node in 
T with two out-edges. The graph G is strongly connected if there is a directed path from 
V to w for any pair of nodes v,w. Any directed path from u to if in a strongly connected 
graph can be extended to a spanning tree rooted at w. Some general references for graph 
theory are [6} and |14j. 
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If G is labeled, then r inherits a labeling from G and we define 

= n 

a 

X — rydr 

Assume that G has no self-loops. We order the node set {vi, . . . ,Vn} of G and let a^j- be 
the label of the edge Vi — t- Vj. Further, we set ajj = for i ^ j ii there is no edge from 
Vi to Vj and Oj^j = 0. Let C{G) = {aij} be the Laplacian of G, that is, the matrix with 
Qjj- = Uj^i if i 7^ j and Oi^i = — Y^2=i ^i,k^ such that the column sums are zero. Any matrix 
whose column sums are zero can be realized as the Laplacian of a directed labeled graph 
with no self-loops. 

For each node Vj, let Q{vj) be the set of spanning trees of G rooted at Vj. Let C{G)(^ij^ 
denote the determinant of the principal minor of C{G) obtained by removing the i-th row 
and the j-th column of C{G). Then, by the Matrix- Tree theorem p?j: 

Note that for notational simplicity we have defined the Laplacian as the transpose of how 
it is usually defined and the Matrix- Tree theorem has been adapted consequently. 

3. Chemical reaction networks 

We introduce the definition of a CRN and some concepts related to CRNs. See for 
instance [9l [11] for extended discussions. 

Definition 3.1. A chemical reaction network (CRN) consists of three finite sets: 

(1) A set 5 of species. 

(2) A set C C M_,^ of complexes. 

(3) A set 7^ C C X C of reactions., such that (y,y) ^ TZ for all y G C, and if y £ C, then 
there exists y' £ C such that either {y,y') £ TZ or {y',y) G TZ. 

Infiow and outfiow of species are accommodated in this setting by incorporating the 
complex E M_,_ and reactions ^ A, A ^ 0, respectively [7j. 

Following the usual convention, an element r = {y,y') E 7^ is denoted by r: y — )• y'. 
For a reaction r: y ^ y', the initial and terminal complexes are denoted by y{r) := y 
and y'{r) := y' , respectively. By definition, any complex is either the initial or terminal 
complex of some reaction. 

Let s be the cardinality of S. We fix an order in S so that S = {Si, . . . , Ss} and identify 
M"^ with M*. The species Si is identified with the z-th canonical vector of with 1 in 
the i-th position and zeroes elsewhere. An element in is then given as ^l^i KSi. In 
particular, a complex y E C is given as y = Yll=iyi^i or (yi, . . . , j/s)- If r is a reaction, 
yi{r),y[{r) denote the i-th entries of y{r),y'{r) respectively. 

Definition 3.2. We say: 

(i) yi is the stoichiometric coefficient of Si in y. 

(ii) If yj 7^ for some i and y £ C, then Si is part of y, y involves Si, and r involves Si 
for any reaction r such that y(r) = y or y'(r) = y' . 
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(iii) Si, Sj G S interact if Ui, yj 7^ 0, i j, for some complex y. 

(iv) y G C reacts to y' £ C if there is a reaction y ^ y' . 

(v) y £ C ultimately reacts to y' £ C (denoted y =^ y') if there exists a sequence of 
reactions y ^ y^ ^ ■ ■ ■ ^ y"^ ^ y' with y^ G C. 

(vi) Si £ S produces Sj E 5 if there exist two complexes y, y' with yi ^ 0, y'^ 7^ and a 
reaction y ^ y' . 

(vii) Si £ S ultimately produces Sj € 5 if there exist Si^, . . . , Si^ with Si^ = Si, Si^. = Sj 
and such that -Sj^..! produces Si^ ior k = 2, . . . ,r. If each S'jj, belongs to a subset 
5q, C 5 for k = 2, . . . ,r — 1, then S'j ultimately produces Sj via 5q. 

If y ultimately reacts to y' then y and y' are linked. Being linked generates an equivalence 
relation and the classes are called linkage classes. Two complexes y, y' are strongly linked 
if both y ^ y' and y' ^ y. Being strongly linked also defines an equivalence relation and 
the classes are called strong linkage classes. 

We introduce an example (which we will refer to as the main example) that we use to 
illustrate the definitions and constructions below. Consider the CRN with set of species 
S = {Si, . . . , Sg} and set of complexes C = {S1 + S2, S4, S1 + S3, S5, 53 + ^4, Sq, S2 + S5, Si + 
Sj, Sj + Ss, Sq, S2 + S3 + Ss}, reacting according to 

-S*! + 5*2 < ^ »S'4 Si + 5*3 < ^ 5*5 S3 + 5*4 < ^ Sq < ^ S2 + 5*5 

I 

5*7 + Ss < ' 5*9 ^ 5*2 + 5*3 + Ss Si + 5*7 

That is, the set of reactions TZ consists of 

ri : Si+ S2 —i- S4 r2 : 54 — > S*! + S2 r^: Si + S3 —?- S5 r^: S5 —s- Si + S3 
rs : 5*3 + 5*4 Sq re : S^ —j- S3 + S4, rj: S2 + S5 —?- Sq rg: Sq—j- S2 + S5 
rg: Sq —J- Si + Sj rio: Sj + Ss —i- Sg rii : Sg —J- Sj + Ss ri2: Sg S2 + S3 + Ss 

This system represents a two substrate enzyme catalysis with unordered substrate binding 
[5] in which Si is an enzyme, 5*2,53 are substrates, 54,55,^6 are intermediate enzyme- 
substrate complexes, and ^7 is considered the product of the reaction system. The product 
dissociates via catalysis by an enzyme Ss and the formation of an intermediate complex Sg. 

The stoichiometric coefficients of all species that are part of a complex are one. The complex 

— 9 

52 + 53 + 58 involves 52, 53, 58 and its vector expression is (0,1,1,0, 0, 0, 0, 1, 0) £ M._^_. The 
species 53 and 54 interact. The complex Si + 53 reacts to the complex 55, implying that 
species Si produces species 55. Also, the complex 53 + 54 ultimately reacts to Si + 57, 
and 57 + 58 ultimately reacts to 52 + 53 + 58. It follows that 53 ultimately produces 57 
and 58- 

4. Mass-action kinetics 

The molar concentration of species Si at time t is denoted by Cj = Ci{t). To any complex 
y we associate a monomial = 131=1 '^T- ^'^^ example, if y = (2, 1, 0, 1) £ M^, then the 
associated monomial is = cfc2C4. 

We assume that each reaction r : y ^ y' has an associated positive rate constant 
ky^yi £ M-i- (also denoted k^.). The set of reactions together with their associated rate 
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constants give rise to a polynomial system of ODEs taken with mass action kinetics: 

(4.1) Ci= ^ ky^yic^iyi-yi), Si£S. 

These ODEs describe the dynamics of the concentrations q in time. The steady states of 
the system are the solutions to a system of polynomial equations in ci , . . . , obtained by 
setting the derivatives of the concentrations to zero: 

(4.2) 0= Yl Wc''(y^-y^), foraUi 
These polynomial equations can be written as: 

(4.3) = krcy^''^y'i{r) - krcy^''^y,{r), for all i. 

It is convenient to treat the rate constants as parameters with unspecified values, that 
is as symbols (as we did in the example in the introduction). For that, let 

Con = {ky^yi\y ^ y £ U} 



be the set of the symbols. Then, the system (4.2) is a system of polynomial equations in 
ci, . . . , Cs with coefficients in the field ]R(Con). 

Only non-negative solutions of the steady state equations are biologically or chemically 
meaningful and we focus on these only. The concept of S-positivity introduced above will 
be key in what follows. Consider the main example and denote by ki the rate constant of 
reaction r^. The mass action ODEs are: 

C'l = -kiCiC2 + /C2C4 - k3CiC3 + k^C^ + kgCQ C5 = fcsCiCs - A:4C5 - A;7C2C5 + ksCQ 

C2 = —kiciC2 + A;2C4 — kjC2C^ + ksCQ + k^CQ cj = kgce — /ciocycs + fcucg 
C3 = -kscics + kiC5 - k^csCi + k^CQ + ki2CQ c's = -ki^c-jc^ + fcncg + ki2Cg 
64, = kiciC2 - k2C4 - k^C2,C4, + keCQ cg = kiocrcg - kucg - k^cg. 

C'e = ^5C3C4 - ksCQ + kjC2C^ - kgCQ - kgc^ 

Take for instance species Si. The only reactions that involve Si are ri, r2, rs, r4, rg. ri,r3 
involve Si in the initial complex and thus the monomials contain ci and have negative 
coefficients. Similarly, r2,r4^,rg involve Si only in the terminal complex and thus the 
monomials do not include ci and have positive coefficients. 

5. Conservation laws and P-semiflows 

The dynamics of a CRN system might preserve quantities that remain constant over 
time. If this is the case, the dynamics takes place in a proper invariant subspace of M*. 
Let X ■ x' denote the Euclidian scalar product of two vectors x, x'. 

Definition 5.1. The stoichiometric subspace of a CRN, {S,C,TZ), is the following subspace 
of W: 

^ = {y - y\y y & 

A semiflow is a non-zero vector uj = (Ai, . . . , A^) G E-*-. If A, > for all i, then a; is a 
P-semiflow. 



8 



Variable elimination in chemical reaction networks 



By the definition of the mass action ODEs, the vector c points along the stoichiometric 
subspace F. The stoichiometric class of a concentration vector c G M_|_ is Cc = {c+r}nM_|_. 
In CRNT, two steady states c, cf are called stoichiometrically compatible if c — c' G F. This 
is equivalent to u; ■ c = uj ■ c' for all u G F"*". 

If w = (Ai, . . . , As) G F"*-, then Yli=i ^i^i — 0- This implies that the linear combination 
of concentrations J2i=i independent of time and thus determined by the initial con- 

centrations of the system. In particular, any steady state solution of the system preserves 
the total initial amounts and lies in a particular coset of F. 

A linear combination Yli=i ^i'^i that is independent of time gives rise to an equation, 
called a conservation law, with a fixed total amount u G M+: 

s 

(5.2) uj = y^^XjCj. 

i=l 

A basis {u^, . . . of F-*- gives a set of independent semiflows and thus a set of in- 
dependent conservation laws: if = Yll=i K^i ^^'^ total amounts oJ^,. . . ,oj'^ G M+ are 
given, we require the steady state solutions to satisfy: cj' = Yli=i ^• 

In the main example, the dimension of F is 5: 

F = {Si + S2 — S4, Si + Ss — S5, S3 + S4 — Sq, Si + S7 — Sq, S7 + Ss — Sg). 
Thus the space F-*- has dimension 4 and a basis is: 

(5.3) cu^ = Si + S4 + S5 + Se co^ = S8 + Sg 

u;^ = S2 + S4 + S6 + Sy + Sg co^ = S2 + S3 + S4 + S5 + 2Se + 2Sy + 2^9. 

The conservation law corresponding to o;^ is cJ^ = ci -|- C4 + C5 -|- cq for a given uJ^ G M+. 

Remark. Questions like "How many steady states does a system possess?" refer to 
the number of steady state solutions that fulfill the conservation laws with the same total 
amounts, or, equivalently, to the number of (stoichiometrically compatible) steady states 
in each stoichiometric class. If this restriction is not imposed and conservation laws exist, 
then the steady state solutions form an algebraic variety of dimension at least 1. 

Remark. Not all CRNs have semiflows. Consider for instance the CRN with s = 6 
and reactions Si ^ S2 , S2 ^ S3, Si + S2 + S3 — )• Sq, S4 + S^ — )■ 5*6, 5*4 — )• S^, and 
S5 Si. The stoichiometric subspace is and thus F-*- = 0. If the last reaction is 
removed, then the stoichiometric subspace has dimension one and there is one P-semiflow: 
UJ = 2Si + 25*2 -I- 2^3 -I- 3S'4 + + QSq. In general, a basis for F-*- consisting of P- 
semifiows is neither guaranteed. Consider the following CRN with s = 3 and one reaction, 
A + B + C ^ A. There is not a basis of F-*- = {A, B — C) consisting of P-semiflows alone. 
This CRN is not biochemically reasonable. 

Lemma 5.4. The following statements are equivalent: 

(i) The stoichiometric class Cc = {c + T} D M^, c G M^, is compact. 

(ii) Fnl+ = {0}. 

(Hi) F-*- has a basis {u^, . . . ,lo'^} of P-semiflows with Aj > if to^ = (A{, . . . , As). 
(iv) There is an element 00 = (Ai, . . . , A^) of with Aj > for all i. 
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Proof. We will prove (i)=^(ii)=^(iii)=^(iv)=^(i). Assume that Cc is compact and consider 
V := rnM^. If y 7^ {O}, then for any non-zero v £ V, the set c + {v) is unbounded in M_|_. 
Hence Cc cannot be compact and thus (ii) must be the case. If (ii), then F-*- n / {0} 
and further T-^ nl+ ^ bd(M+). If the latter was not the case, then also rnbd(M+) / {0}, 
contradicting (ii). Hence, there exists an open set C n in F"*", and we can choose 

a basis {(^^, ■ ■ ■ , w"^} of P-semiflows with lo^ S O, that is, > 0. Thus (iii) is fulfilled, (iii) 
gives (iv) directly. Assume (iv). For x = {xi)i G Cc, UJ ■ x = u ■ c is independent of x. Since 
Aj > and Xi > 0, Xi < {u ■ c)/\i for all i and thus Cc is bounded. Since it is a closed set, 
Cc is compact and (i) is proven. □ 



Lemma 5.4 is well-known in dynamical systems theory and Petri Net theory. In the 
latter semiflows are known as P-invariants (place invariants) j21j . 

Remark. All conservation laws might not be obtained from semiflows [12], that is, 
the semiflows in F"*- might not give the minimal afiine space in which the dynamics of 
the system takes place. There can be additional conservation laws depending on the rate 
constants and not merely on the stoichiometric coefficients. The next lemma is proven in 
[12\ and stated here for future reference: 

Lemma 5.5 ([12j, §6). // each linkage class contains exactly one terminal strong linkage 
class, then all conservation laws correspond to semiflows. 

As shown in |i2j , any weakly reversible network fulfills the condition of the lemma. Also, 
the main example fulfills the criterion. The CRN with reactions ri : 5*1 — )• 5*2 , r2 : Si — )■ S3 
and r^: S2 + Ss ^ 2Si does not fulfill it [3]. Here, F"*- = {Si -|- 5*2 + S3) providing the 
conservation law ci + C2 + C3 = uj. However, when ki = k2 = k^, then ci -|- 2c2 is also 
conserved. 

Minimal and terminal semiflows. 

Definition 5.6. The support of a semiflow u = (Ai, . . . , A^) is the set S{uj) = {Si\ Aj 7^ 0}. 
We say that w is 

(i) minimal if for any semiflow uj with S{oj) C S{uj), there is a G R such that auj = ui. 

(ii) terminal if any semiflow u with S{ljj) C S{u}) satisfies 5(a3) = S{uj). 

That is, a semiflow w is minimal if any semiflow given by a linear combination of the 
species in its support is a multiple of u and terminal if there is no semiflow with smaller 
support. 

Lemma 5.7. (i) A semiflow is minimal if and only if it is terminal, 
(ii) If uj is a P-semiflow that is not minimal, then there is a P-semiflow uj such that 
S(uj) C S{uj). 

Proof, (i) If is a minimal semiflow then by definition any semiflow cu with S{oj) C S{uj) 
satisfies oj = auj for some a G M. Thus, S{cjj) = S{u}), which implies that uj is terminal. To 
prove the reverse, assume that oj is terminal but not minimal, that is, there exists cj such 
that S(u}) = S{lj) and tj 7^ ao; for all a G M. Let I = {i\Si G S(uj)}, oj = (Ai, . . . , As), and 
OJ = (Ai, . . . , Xs). Choose u £ I such that |Au/Am| > |Ai/Aj| for all i £ Z and define 7 = A„ 
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and 7 = Au. Then 

s s 

Q := -yuj - juj = ^(A„Ai - Xu\i)Si = ^ fiiSi 

i=l i=l 

is a semiflow, since tD 7^ (otherwise lj = auj for some a). Since fiu = 0, S{uj) C S{u}), 
which contradicts that lv is terminal. 

(ii) If is a P-semiflow that is not minimal, then there exists a semiflow u such that 
S{u}) C S{ljj). The construction above provides a new semiflow lD. Since Aj > for all 
I G X, we have lA^IAj — Au|Aj| > and either A^ > and /ij > for all i, or A^j < and 

< for all i. Hence, either u) or —u) is a P-semiflow fulfilling (ii). □ 

Therefore, there cannot exist two linearly independent minimal P-semiflows with the 
same support. For example, if S'i+S'2+5'3 is conserved and minimal, then AiS'i+A25'2+A3S'3 
with Ai 7^ A2 cannot be a semiflow. We will see below that the P-semiflows uj^,oj'^,uj^ in 



(5.3) of the main example are minimal. However, co^ is not minimal since S{uj^) C S{u^). 

The species graph does not characterize the CRN uniquely, since information coming 
from the stoichiometric coefficients is ignored. For instance, the following two systems have 
the same species graph: 

(5.8) Ri = {A + B ^2C,C ^ A}, R2 = {A + B ^ C,C ^ A}. 

6. Species graph 

Given a CRN {S,C,TZ), we define the species graph Gs as the labeled directed graph 
with node set S and a directed edge from Si to Sj with label r: y ^ y' whenever yi 
and y'j ^ 0. That is, there is a directed edge from Si to Sj if and only if Si produces Sj. 
There can be multiple edges with different labels between a pair of nodes. In addition, if Si 
is involved in the initial and the terminal complexes of a reaction, then there is a self-edge 
Si ^ Si. The species graph of the main example is depicted in Figure [l} 

Remark. A reaction r: y — )• y' is called reversible if the reaction y' ^ y also exists. 
In the main example, all reactions but rg,ri2 are reversible. In contrast to other papers 
[Tl[23], we consider reversible reactions as two (independent) irreversible reactions. Thus, 
reversible reactions provide two edges with opposite directions and different labels in the 
species graph. This is required when we consider spanning trees in Section [8] Changing a 
reaction from being reversible to irreversible does not change the stoichiometric subspace 
and a system with all reactions considered irreversible has the same (P-)semiflows as a 
system with some (all) reactions considered reversible. However, the steady states might 
depend on whether reactions are reversible or not. 

Definition 6.1. A graph G with node set Sa is a subgraph of Gg if 5q, C 5 and the labeled 
directed edges of G are inherited from Gg. We denote G = Gs^. Further, 

(i) Sa is full if any reaction involving Si G Sa appears at least once as a label of an edge 
in Gsa- If this is the case, then is said to be full. 

(ii) Sa is non-interacting if it contains no pair of interacting species and all stoichiometric 
coefficients are either or 1, that is, = 0, 1 for all Si E Sa and y G C. If this is the 
case, then Gg^ is said to be non-interacting. 

(iii) If Sa is full and non- interacting, then Sa is a cut of S. 
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Figure 1. Species graph of the main example. 

The definition of subgraphs of Gs extends to subgraphs of . We depict in Figure [2] 
four different subgraphs of the species graph of the main example, corresponding to four 
different subsets 5q, C 5. 

The proof of the following lemma is left to the reader. 

Lemma 6.2. Let 81,82 ^ 8. Then 

(i) If Gsi and are full, then so is GsiuS2- 

(ii) If GsiVjS2 'i'^ non-interacting, then so are Gs^ and Gs^- 

If Gsi n Gs2 = and GsiuS2 = ^^^i U Gs2, then the reverse statements are also true. 

It follows that if 8a is a cut, then the node set of any connected component of G^^ is 
also a cut (as illustrated in Figure [2]||a) ) . The next lemma connects some properties of full 
and non-interacting graphs that will be used in the sequel. A non-empty subgraph G^^ of 
Gs^ is proper if Gs^ + Gs^- 

Lemma 6.3. Let 8a 8 be a subset. 

(i) If Gsa is non-interacting, then any reaction label r appears at most once in Gs^. 
(ii) If Gsa is non-interacting and Si £ 8a is involved in a reaction r that is a label of an 

edge ofGs^, then the edge is to/from Si. 
(Hi) If Gsc is full and Si G 8a is involved in a reaction r, then there is an edge to/from 
Si labeled r. 

(iv) If Gsa is full, has no repeated reaction labels and the stoichiometric coefficients of its 
nodes in all complexes are either or 1, then 8a is a cut. 

(v) If Gsa has no repeated reaction labels and Gs^ is connected, then Gs^ has no proper 
full subgraphs. 

Proof, (i) Assume that a reaction r appears in two different edges Si A Sj and Su Sy 
of Gsa ■ li Si ^ Su or Sj / Sy , then either Si and Su or Sj and Sy interact and thus Gs^ 
is not non-interacting. 

(ii) -(iii) In both cases there is an edge in Gs^ labeled r: In (ii) by assumption and in 
(iii) because Gs^ is full. Assume that the edge with label r is not from/to Si but between 
Sj, Su £ 8a, j, u ^ i. Then, Si and Sj (or Su) interact reaching a contradiction in case (ii) 
and implying that there is an edge with label r between Si and Su (or Sj) in case (iii). 
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(a) cS^ = {81,84,, 85,8(1, 58, 59} 




^6 > Si H '± Sg 



(c) Si = {S4, 55,56,57,59} 




(b) 5q = {52, 54, 56, 57, 59} 




(d) St = {82,83,84,85,86,87,89} 



Figure 2. Species subgraphs of the main example. The node sets in (a) 
and (b) are cuts but the graph in (a) has two connected components and 
in (b) one; (c) the node set is non-interacting but not full and cannot be 
extended to a cut; (d) is a full graph but two of its nodes, S2, S3 interact. 



(iv) Assume that there are two (different) interacting nodes Si,Sj G Sa- Then there 
exists a reaction r: y ^ y' such that yi = yj = 1 or y'- = y'j = 1. Let us assume that 
Vi ~ Vj ~ ^ ^'^^ other case follows by symmetry. Since Gs^ is full, r is the label of an 
edge in Gs^- Let Su (potentially equal to Si or Sj) be the end node of the edge. Since 
S'j, Sj, Su ^ Sa, the edges Si A Su and Sj A Su are in G^^, contradicting that there are 
no repeated labels. Hence, Sa is non-interacting and hence a cut. 

(v) Assume that there is a proper subgraph G. Since Gs^ is connected one can find a 
node Si in Gg^ that is not in G, and a node Sj in G for which there is an edge Si A Sj or 
Sj Si in Gsa ■ By assumption a label appears at most once in Gs^ ■ Thus G is not full 
since r is not a label of an edge in G, but r involves Sj £ G. □ 



It follows from Lemma 6.3' i,iii) that if Sa is a cut, then all reactions involving Si G Sa 
are edges of G^^ to/from Si and they appear exactly once. From (i,v) we find that non- 
interacting connected graphs have no proper full subgraphs. In particular, if Sa is a cut 
such that Gsa is connected, then Gs^ has no proper full subgraphs. 



7. Semiflows and the species graph 

In this section we explore the relationship between (P-)semiflows to full subgraphs of 
Gs- We come to the main results on semiflows in relation to variable elimination: (1) any 
semifiow whose support is a cut and the associated graph is connected has all non-zero 
coefficients equal; (2) the support of a semifiow is non-interacting if and only if it is a cut. 
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Lemma 7.1. Let u be a semiflow. 

(i) If uj is minimal, then Gg(^^'j is connected. 

(a) If UJ is a P-semiflow or S{uj) is non-interacting, then G^^^-^ is full. 
(Hi) If Lo is a P-semiflow or S{uj) is non-interacting, and G^^^^^ has no proper full sub- 
graphs, then UJ is minimal. 

Proof, (i) The idea is that if Gs(^^-j is not connected, then the semiflow splits into two, 
contradicting minimality. If is not connected, then Gg(^^-j = Gi U G2 with Gi,G2 

being two non-empty disjoint subgraphs of Gs(^^^y Let Si, S2 be the node sets of Gi,G2, 
respectively. If cj = Y2i=i ^i^i, let Uk = Z^i|5-g5j. ^iSi 7^ 0, fe = 1, 2. Since 5i n52 = 0, w = 
UJI+UJ2- By hypothesis, = for all y — )• y' e 7^ and thus = wi •(?/'— y)+a;2-(y'—y). 

Since Gi, G2 are disjoint, there is no edge between a node in Si and a node in ^2. If y — )• y' 
is a reaction with y'^ — yi^O for some Si G Si, then yj = y'j = for all Sj E ^2 and trivially 
0J2-{y'—y) = 0. Thus a;i-(y'— y) = 0. By symmetry, we have that a;i-(y'—y) = a;2-(y'— y) = 
for all reactions y ^ y' and hence uok is a semiflow for k = 1,2. 

(ii) Let Si € S{uj) and r: y — )• y' be a reaction with either yj or y[ / 0. We assume that 
yi and the other case follows by symmetry. We want to show that r is a label in G^^^), 
that is, y'j / for some Sj G S{uj). By hypothesis, uj • {y' — y) = 0. If w = (Ai, . . . , A^) is a 
P-semiflow (that is, Aj > 0), we have uj ■ y' = uj ■ y > Xiyi > 0. If S{uj) is non-interacting, 
then yj = for any i ^ j such that Sj G S{uj) and hence uj ■ y' = uj ■ y = Xiyi / 0. Either 
case Xjy'j 7^ for some j and hence y'j ^ for some Sj £ S{uj). Note that the case j = i is 
accepted. 

(iii) Assume that uj is not minimal. Then by Lemma |5.7| ^i) there exists a semiflow uj 
such that S{uj) C S{uj). If w is a P-semiflow, then by Lemma |5.7[ ii) we can assume that 
w is a P-semiflow. If S{uj) is non- interacting, then so is S(uj). Using (ii), Gs(^^) is full and 
thus a proper full subgraph of reaching a contradiction. □ 



Consider the P-semiflows a;* (5.3) of the main example and the subsets 5^ in Figure 
21 Here, 5(0;^) U S{uj^) = S^, Sju?) = S^ and S{uj'^) = S^, giving fuh subgraphs. The 



7.1 



iii) that uj^,uj'^,uj^ are 



two components of G51 and the graph G52 have no proper full subgraphs, while G52 is 
a proper full subgraph of Thus, it follows from Lemma 

minimal. 

Lemma [7.1[ iii) cannot be reversed. Consider for example the CRN R2 in (5.8). The vec- 

ri n 

tor UJ = A-\-B-\-C is a minimal P-semiflow. The associated graph G^(^^-j is A < _ ^ C B 

ri 



>"2 



and has a proper full subgraph A 



C ■ Further, W2 = ^ + C is a P-semiflow for the 



CRN Ri in (5.8), but there is not a P-semiflow involving all species. This implies that if 
G is a full connected subgraph of G and G corresponds to a P-semiflow, then G does not 
necessarily correspond to a P-semiflow. Similarly, if G corresponds to a P-semiflow, then 
G does not necessarily correspond to one too. 



Proposition 7.2. Let Sa S be a subset. If Sa is a cut, then cj = ^ 
P-semiflow. In this case, uj is minimal if and only if Gs^ is connected. 



Si is a 
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Proof. The stoichiometric coefficients of Si G Sa are either or 1 since Sa is a cut. Consider 
r : y — )• y' E 7^. If ?/j , y • = for all i such that Si G 5^, then clearly u)-{y' — y) = 0. Since 5q, 
is non-interacting, if 7^ for some Si G Sa, then yj = for all Sj £ Sa such that i ^ j. 



From Lemma 6.3 }i,iii), r is the label of exactly one edge of Gs^ connected to Si, and there 
exists exactly one species Sj G Sa such that y'j ^ 0. Thus, uj- {y' — y) = yi — y'j = 1 — 1 = 0. 
This shows that u; is a P-semiflow. 

For the second part of the statement, Lemma |7.1[ i) implies that if uj is minimal then 



Gs^ is connected. On the other hand, if Gs^ is connected, then by Lemma 6.3 l,v), 
has no proper full subgraphs and by Lemma [7.1[ iii), oj is minimal. □ 

The reverse is not true: In the reaction system Si + S2 ^ Yi + Y2, S3 ^ Si, Yi ^ S3, 
Y2 ^ S3, the vector w = 5i + ^3 + Yi + 12 is a P-semiflow but the set 53, 1^, 12} is 
not a cut. Further, to is minimal. Thus, the non-interacting property cannot be read from 
the coefficients of the P-semiflow. 



Remark. If Sa is a cut, we denote the P-semiflow given in Proposition 7.2 by uj{Sa) = 
Si|5ie-SQ '^'^^ operation uj{-) defines a map between the set of cuts and the set of 
P-semiflows such that sets with connected associated graphs are mapped to minimal P- 
semiflows. The operation S{-) defines a map between the set of P-semiflows and the subsets 
of S with full associated graphs such that minimal P-semiflows are mapped to subsets with 
connected associated graphs. The map S{uj{-)) is the identity. 

From Lemma |7.1[ ii) and Proposition |7.2| we derive the following corollary. 

Corollary 7.3. (i) Let lo be a semiflow such that S{uj) is non-interacting. Then S{uj) 
is a cut. 

(a) If Sa S is a non-interacting subset and uj{Sa) is not a P-semiflow, then Sa is not 
a cut and there is no semiflow with support Sa. 

Therefore, there is a one-to-one correspondence between cuts and P-semiflows with non- 
interacting support and non-zero entries equal to one. The results of this section give rise 
to the following corollary. 

Corollary 7.4. Let Sa be a non-interacting set such that Gs^ is connected. 

(i) If Sa is a cut, then any semiflow with support included in Sa is a multiple of uj{Sa). 
(ii) If Sa is not a cut, then there are no semiflows with support included in Sa- 



Proof, (i) Follows from Proposition 7.2 and Lemma 5.7 (ii) Follows from Lemma 6.3 l,v), 
Lemma |7.1[ ii) and the fact that Gs^ is not full. □ 

Consider the non-interacting subset 5^ of the main example. The vector to = Si + S5 + 
Sq + Sj + Sg is not a P-semiflow since u ■ {Si + S2 — S4) = — 1 / 0. It follows that 5^ is 
not a cut and there is no semiflow with support included in 5^. Reciprocally, consider the 
P-semiflow lo^. Since Sa = S{uj^) is non-interacting, it is a cut. Further, since the graph 
G52 is connected and 5^ is a cut, any semiflow involving the species in 5^ is a multiple 
of w^. Checking if a set is non- interacting is straightforward from the set of complexes. 
However, checking that the associated subgraph is full can be tedious. We have shown 
that the relationship between semiflows and cuts gives easy conditions for determining if 
a non- inter acting set is a cut. 
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Remark. Using Lemma 6.2 we find that if Sa is a cut but Gs^ is not connected, then 
Sa decomposes into a disjoint union of cuts, Sa = S^U ■ ■ - USJ^, such that G^i^ is connected 
for all i. Thus, a; (5a) = w(5^) + • • • + w(5^) and the P-semiflow uj{Sa) decomposes into 
a sum of minimal P-semiflows. Further, if Gs^ is not connected, e.g. has two connected 
components G^i and G52 , then a; (5^) — a;(5^) is a semiflow with support in Sa = 5^ U5^, 
but it is not a multiple of Lj{Sa)- It follows that requiring Gs^ to be connected is a necessary 
condition in Corollary |7.4[ 



8. Elimination of variables 
Let 5a C 5 be a non- interacting subset. Let Gs^ be the species graph associated to 5q 



and assume that it is connected. By Corollary 7.4, either Sa is a cut and Lo{Sa) is minimal, 
or Sa is not a cut and there is no semiflow with support included in Sa- In what follows we 
discuss conditions such that the concentrations of the species in Sa can be fully eliminated 
from the steady state equations. The discussion depends on whether Sa is a cut or not. 

For simplicity, we assume that Sa = {Si, . . . , Sm}- Let TZa be the set of reactions not 
appearing as labels in but involving some Si £ Sa, and Tla outi''') '^aini''') sets 
of reactions in TZa involving Si G Sa in the initial and terminal complexes, respectively. 
Clearly, TZ'^ = IJI^i ^^a.outi^) U m(^)- Note that VJa = if and only if Sa is a cut. 



We restrict our attention to the steady state equation (4.2) for a fixed Si S Sa- Using 



the expression in (4.3) and the fact that the stoichiometric coefficients are one, we can 



write this equation as Si = Xi — Yi with: 

m m 

-Xi — ^ ^ ^ ^ kj-c^^ ^ ~t~ ^ ^ k'[-c^^ ^ , Yi — ^ ^ ^ ^ krfC^^ ^ ~t~ ^ ^ k^c^^ ^ , 

where the first summand in each term is taken over the edges in Gs^ . Recall that there can 
be more than one edge between two species. Further, since the stoichiometric coefficient 



of Si is or 1, any edge Si — )■ Si provides no summand in equation (4.2). 

Let Ca = C{Sa) = {ci\Si G Sa} = {ci, . - - ,Cm}- Each of the monomials c^^^^ in Yi 
involves Cj, and if another Cfc is involved, then Sk interacts with Si (and in particular 
k ^ {1, . . . , m}). Similarly, c^^^^ in Xi involves the variable Ck if and only if Sk produces 
Si- Further, for r £ T^aini''')' '^^^''^ does not involve any Ck G Ca, while for Sj A Si with 
Sj £ Sa, the only such variable is Cj. It follows that the system is linear in Ca with 
coefficients in R[ConU C"^(5q)], where 

(8.1) C^(5q,) = {ci\Si £ S \ Sa interacts with or produces some Sj G Sa}- 

We write = C^{Sa) for short and note that Cq n = 0. Let yj{r) denote the vector 
in M''"^ obtained from y{r) by removing the j-th coordinate. We have shown that Xi,Yi 



can be written so that equation (4.2) for Si E Sa becomes 



^•2) ^ = X] + ^« 
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where an = ei + di and 



Let A = {oij} be the m x m matrix with aij defined as above, d = {di, . . . , dm) and 
z = (zi, . . . , Zm)- Note that Sa is a cut if and only if z = d = 0. The discussion above 
provides a proof of the following lemma: 



Lemma 8.3. The steady state equations (4.2) for Si £ Sa form an mxm linear system of 



equations in Ca, Ax + z = 0, where the entries of the matrix A and the independent term z 
are either zero or S-positive in M[ConU C^]. Further, Sa is a cut if and only if z = d = 0, 
in which case the system is homogeneous. 

If A has maximal rank m, then the system has a unique solution in M(ConUC^). By 



Corollary 7.4 , if Sa is a cut then the column sums of A are all zero and the system cannot 
have maximal rank. If Sa is not a cut then there are no semiflows with support in Sq^. The 
column sums of the matrix A are (for column i): 

m m 

^ak,i = Y, E Kcy'^'^-Yl E krcy^^'^ + di = d,. 

k=l i^k i = i s-^s- 

These are zero as polynomials in R[ConUC^] if and only if 'R-aouti''') — ® for all and 
the condition d = is equivalent to the column sums being zero. If Sa is not a cut, but 
c? = 0, then z 7^ as a tuple with entries in ]R[ConUC^]. It follows that the system is 
incompatible in M(ConUC^), because -Zj 7^ 0. The only possible non-negative steady 
state solutions must satisfy = for some i and hence Cj = for some Cj G such that 
Sj produces Si £ Sa- 

We proceed now to discuss the case in which Sa is a cut and the case in which it is not 
a cut. Both cases could be merged into a single approach, but the discussion of the first 
situation becomes more transparent when it is treated separately. 

Elimination of variables in a cut. Let 5q C 5 be a cut such that Gs^ is connected. 



For Si G Sa the equations in (4.2) form an m x m homogeneous linear system of equations 



with variables Ca and coefficients in M[ConUC^]. Using equation (8.2), equation (4.2) 
becomes 

m 

(8.4) = ^OjjCj, i = l,...,m. 

Because the column sums of A are zero, A is the Laplacian of a labeled directed graph 
with node set 5^ and a labeled edge Sj — ^ Si, whenever Oij / 0, i 7^ j. Note that 
Oij £ M[ConUC^] is S-positive. We have that Gs^ is (strongly) connected if and only if 
Gsa is. The two graphs differ in the labels and in that multiple directed edges from Si to 
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Sj in are collapsed to a single directed edge in Gs^- Further, the graph Gs^ has no 
self- loops (that is, those of Gs^ are removed by construction). 

By the Matrix- Tree theorem, the principal minors ^(ij) of vl = C{Gs^) are 

Am) = (-ir^'^^""' E ^(^)- 

re0(Sj) 

Thus, A has rank m — 1 if and only if there exists at least one spanning tree in Gs^ rooted 
at some Sj, j = 1, . . . ,m. The next proposition follows from the discussion above. In 
particular, it holds if Gs^ (or equivalently Gs^) is strongly connected. 

Proposition 8.5. Assume that Sa is a cut such that Gs^ is connected and let uJ = X^I^i 
be the conservation law obtained from the P-semiflow uj{Sa)- The following statements are 
equivalent: 

(i) UJ = X^I^i ^-^ ^''^^y conservation law with variables in Ca- 
(ii) Gsa ^^-5 fli least one rooted spanning tree. 
(Hi) The rank of A is m — \. 

Since is connected and Sa is a cut, any semiflow with support in Sa is a multiple of 
uj{Sa)- The proposition says that Gs^ has a rooted spanning tree if and only if there are 
no other conservation laws with concentrations only in Ca- 

Remark. Let Cq be the set of complexes involving at least one species in Sa. Consider the 



linkage classes in Ca given by the relation "ultimately reacts to" (Definition 3.2). If there 
is a rooted spanning tree, then the root must be in a terminal strong linkage class, because 
the elements in such a class cannot react to complexes outside the class. Further, using the 
same reasoning, there cannot be two terminal strong linkage classes. Consequently, if there 
is a rooted spanning tree, there is only one terminal strong linkage class. This remark is 



closely related to Lemma 5.5 



For simplicity we assume that there exists a spanning tree rooted at Si. Then, the 
variables C2, . . . , can be solved in the coefficient field M(ConU U {ci}). In particular, 
using Cramer's rule and the Matrix- Tree theorem, we obtain 

(8.6) c, = = ^4§T^i = V'^^^")^!' ^^■(^") = E ^(^) 

and j = 1, . . . ,m. Since there is a spanning tree rooted at Si, it follows that ai{C^) is S- 
positive and o-j{G^) is either zero or S-positive in ]R[ConU C^]. If the graph Gs^ is strongly 
connected, then (Tj{C%) ^ for all j and any choice of root Sj could be used instead of 
Si- The arguments given above and the definition of aij provide a proof of the following 
lemma. 

Lemma 8.7. If Ck G is a variable of the function crj(C^) for some j, then there exists 
Si G Sa that interacts with Sk and Si ultimately produces Sj via Sa. Specifically, there is 
a complex involving Si and Sk, and a complex y^ involving some species Su ^ Sa, such 
that y^ reacts to and Su ultimately produces Sj via Sa- If Gs^ is strongly connected, 
then the reverse is true. 
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The sum of the concentrations in Ca is conserved. Using the equation a; 
obtain 

a; = (l + (^2(C^) + --- + '^m(C^))ci, 
where the coefficient of ci is S-positive in M(ConUC^). Thus, 

uj uj(Ji{C%) 



Em 
i=iQ, we 



with Tp-^ being an S-positive rational function in with coefficients in M(ConU{Zi;}). 
Observe that a; becomes an extra parameter and can be treated as a symbol as well. 
Further, if w is assigned a positive value, then ci > at steady state for positive values of 
Cq. By substitution of c\ by obtain 

(8.8) Cj = (^,(C7^) := <^,(C^)^i(C^), j = 2, . . . , m, 

with Tp- being either zero or an S-positive rational function in (7^ with coefficients in 
M(ConU{tJ}). 

Proposition 8.9. Let Sa Q S be a cut such that Gs^ is connected. Assume that there is 
a spanning tree of Gs^ rooted at some species Si . Then, there exists a zero or S-positive 
rational function ipj in with coefficients in M(Con), such that equation (4.2) for Cj E Ca 
is satisfied in R(ConU C^) if and only if 



C,;, 



Cj G Ca • 



Further, there exists an S-positive rational function ip^ in C^ with coefficients mM(ConU{a;}), 
such that the conservation law uJ = X^^^ is fulfilled if and only if Ci = ^j(C^). 



Consider the main example and the cut Sa = {Si, S^, S5, Sq} corresponding to a con- 
nected component of the graph C^i in Figure [2Fa). System (8.4) becomes 



v 



-kiC2 - hcs 
kiC2 





k2 

-k2 - k^cz 


^5C3 



^6 

-ke - ks- kg J 

The column sums are zero because of the conservation law uj^ 



ki 


A;4 - kYC2 
/C7C2 



0. 



C4 

V C6 y 

_ ci -I- C4 + C5 -I- C6. The 

graph Gsa is strongly connected (Figure |3|^ a)), as is observed in many real (bio)chemical 
systems. Thus rooted spanning trees exist and the system has rank 3. This also follows 
from Proposition 8.5 and Lemma 5.5, since each linkage class of the CRN has exactly one 
terminal strong linkage class. 
The polynomials CTj are: 

fji = k2k4{kQ + /cs + kg) -\- k2kj{kQ -\- kg)c2 + k/ik^{k^ + kg)c^ + k^kjkgC2Cz 

(T4 = kikiik^ + ^8 + kg)c2 + kikj{kQ + kg)c2 + k^kQk-jC2Cz 

0-5 = k2kz{k& + ks + kg)c3 + ksk^iks + kg)cl + kik^ksC2Cz 

(T6 = {kik4k5 + k2k^kj)c2Cz + kik^kjc^c^ -|- k^k5k7C2c\. 

Each monomial in Uj corresponds to a spanning tree rooted at Sj . The species 5*2 , ^3 
are the only species interacting with a species in Sa and thus only C2 , C3 appear in the 
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5, ^ 



(a) Sa = {5i, 54, 55, Se} 




(b) 52 = {54,55,56,57,59} 



Figure 3. (a) The graph Gs^ in the main example for the cut Sa and (b) 
The graph for the non-interacting set S^, which is not a cut. 



Using (8.6) and (8.8) we find the steady state expressions of ci, 04,05,05 in 
the total amount uJ^, and the concentrations 02,03. 

by computing the principal minors of A 



expressions. 

terms of the rate constants, 

Remark. It is straightforward to find aj 
using any computer algebra software. The advantage of the Matrix- Tree description in the 
theoretical discussion is that S-positivity of the solutions is easily obtained. 



Elimination of variables in a subset that is not a cut. Let 5^ C 5 be a non- 
interacting subset that is not a cut and assume that Gs^ is connected. As discussed above, 
if the column sums are zero then there are no positive steady state solutions. If the column 
sums are not all zero, then the matrix A is not a Laplacian. However, A can be extended 
such that its determinant is a principal minor of a Laplacian. 

Consider the labeled directed graph Gs^ with node set 5^ U {*}. We order the nodes 
such that Si is the i-th node and * the (m + l)-th node. The graph Gs^ has the following 
labeled directed edges: Sj — '-^ Si if ajj / and i / j. Si — ^ * if / 0, and * ^ 5j if 
Zi 7^ 0. All labels are S-positive in R[ConU C^]- Let £ = {Xij} be the Laplacian of Gs^. 



li i,j < m, then A 



j. The entries of the last row are Xm+i,i = —di for i < m and the 



entries of the last column are Aj^m+i = Zi for i < m. We conclude that the (m + 1, m + 1) 
principal minor of C is exactly A and thus, by the Matrix- Tree theorem, we have 

(m+l,m+l) = X] ^('^)- 

Tee(*) 



-i)™det(^) = i-irc 



If there exists at least one spanning tree rooted at *, then (— l)™'det(^) is S-positive in 
M[ConUC^]. In this case the system Ax + z = has a unique solution in M(ConUC^). A 
spanning tree rooted at * exists if and only if for all species Si G Sa, there exists a reaction 
y ^ y' such that y' does not involve any species in 5q,, y involves some Su G Sa and Si 
ultimately produces Su. The existence of such a spanning tree ensures that 71% outif) 7^ ^ 
and thus di for some i. 

Since A is non-interacting but not a cut, there are no semifiows with support in Sa, 
Corollary |7.4rii). Similarly to Proposition 8.5 we obtain the following proposition: 



Proposition 8.10. Assume that Sa is a non-interacting set that is not cut and such that 
Gs^ is connected. Then A has maximal rank if and only if there exists a spanning tree 
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rooted at *. Further, if A has maximal rank then there are no conservation laws in the 
concentrations in Ca- 

If A does not have maximal rank, then there is a vanishing linear combination of the rows 
of ^, = X^fcLi ^kO^kj for all j. If there are no conservation laws in the concentrations in 
Ca, then / J2k=i ^kCk = EfcLi EJli ^k{ak,jCj + Zk) and it follows that YJ^=i ^kZk + 0. 



If Afc > for all fc, then we conclude that the system (8.2) is incompatible in M(ConU C^) 
and there are no positive steady states. 

Assume that a spanning tree rooted at * exists. For z = 1, . . . , m, let cjj be the following 
polynomial in C^, 

re0(S,) 

which is either zero or S-positive in M[ConUC^]. By Cramer's rule, we have 



(~l)™'^(m+l,m+l) ^i^a) ' 

which is either zero or S-positive in M(ConUC^). If there exists at least one spanning 
tree rooted at Si, then cTj 7^ as a polynomial in M[ConUC^]. A necessary condition for 
CTj 7^ is the existence of a directed path from * to Si, which implies that Si is ultimately 
produced from some species 5^ £ S \ Sa- In particular, if Gs^ is strongly connected then 
all concentrations are non-zero as elements in M(ConUC^). 

Consider the set Sa = = {S4, S5, Sq, S^, Sg} in the main example. It is non- 
interacting, not a cut, and G^^ is connected. Further, all species ultimately produce 
Sg E Sa and 5*9 reacts to 52 + + Sg which does not involve species in Sa- Hence a 
spanning tree rooted at * exists. The graph Gs^ is depicted in Figure js^b) and is strongly 
connected. We have that z = {kiCiC2, k^cics, 0,0,0), C^ = {ci, C2, C3, cs} and 

a =kioki2{k2k4,{kQ + kg + kg) + k2kQkjC2+ 

k^k^kgcs + kgikik^cs + k2k7C2 + k^ok7C2C3))cs 

0-4 =kwki2{kiki{kQ + ks + kg) + kikjike + kg)c2 + k3kQk7C3)ciC2Cs 

0-5 =kioki2{k2k3{kQ + ks + kg) + kikr^ksC2 + k^k^^k^ + kg)c3)ciC3C8 

=kiQki2{kikik^ + k2k3k-j + kik^kjC2 + k3k^kjC3)ciC2C3Cs 

0-7 =A;9(A;ii + ki2){kikik^ + k2k3k'j + kik^k-jC2 + k3k^kjc-i)ciC2C-i 

ag =kgkioikik4k5 + k2k3k7 + kikr^kjC2 + k3kr^kjOi)ciC20iC^. 

The concentration ci is only in the label of out-edges from * and thus ci is not in a. 

Proposition 8.11. Assume that there is a spanning tree of Gs^ rooted at *. Then, there 
exists a zero or S-positive rational function ipi in G^ with coefficients in ]R(Con), such that 



equation (4.2) for Ci S Cq, is satisfied in M(ConUC^) if and only if Ci = (^j(C^). 



Remark. The procedure outlined here can be stated in full generality: Consider a square 
linear system of equations Ax + z = 0, such that the entries of z and the off-diagonal entries 
of A are positive and the column sums of A are zero or negative. If A has maximal rank, 
then the unique solution of the system is non- negative. 
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Remark. If the matrix A does not have maximal rank, then we can always selecxt 
a subset of Sa such that the corresponding matrix has maximal rank and proceed with 
elimination of the variables in the subset. 

Remarks. Let 5^ C 5 be any non-interacting subset such that Gs^ is connected. We 
have proven that for all Ci £ Ca, there exists a rational function ipi such that Cj = (pi{C%) 
at steady state, provided some spanning trees exist. When Sa is a cut, the P-semiflow 
oj{Sa) is required in the elimination, while when Sa is not a cut, variables are eliminated 
using only the steady state equations. 

If is not connected, then the results above apply to the connected components 
separately, since the underlying node set of each connected component is non-interacting. 
Further, let 5^,5^ be two non-interacting sets such that Ggi and G52 are disjoint and 
connected. One easily sees that G{Sl)<r^G''{Sl) = %, that is, both sets of variables G{Sl) 
and C(5^) can be simultaneously eliminated. Additionally, if we let Sa = Slvj Si then 
C^{Sa) = C'^{S^) U C"^(5^). For instance, consider Sa = {Si, S^, S5, Sq, Ss, Sg} in Figure 
^a). The associated graph has two connected components, which are strongly connected. 
The concentrations Cj S Ca can be expressed as S-positive rational functions in = 
{c2, 03,07}. 

The conditions to apply the variable elimination procedure are summarized in Table [T] 
The procedure guarantees that if positive values are assigned to all Cj £ C^, then q is 
non-negative. For Cj to be positive, that is, ai ^ 0, the existence of at least one in-edge to 
Si is necessary. Otherwise the concentration at steady state of Si is zero, which is expected 
if Si is only consumed and never produced. Further: 

Proposition 8.12. Let Sa be a non-interacting subset of S such that the concentrations 
Ca can be eliminated from the steady state equations. Each component of the graph Gs^ 
is strongly connected if and only if any steady state solution satisfies cj > for all Cj £ 
Ca ( and for any positive total amount if appropriate ), whenever the variables in C% take 
positive values. 

9. Steady state equations 

Let Sa be any non-interacting subset such that Cs^ is connected and that the variables 
in Ca can be eliminated by the procedure above. Let $u(C^) = be the equation ob- 
tained from Cu = 0, u = m + 1, . . . , s, after elimination of variables in Ca and removal 





Sa cut 


Sa not a cut 


Characterization 

Elimination of 
Ca works if in 
Gs^ ••• 

Table 1. Sun 


00 {Sa) semiflow or z = d = 

3 rooted spanning tree (equiv- 
alent to cJ = Ylc^i^Ca being 
the "only" conservation law in 

Ca) 

imary of the conditions requirec 


J semiflow or z — d 7^ 

3 spanning tree rooted at * 
(implies ^ conservation law in 
Ca and d 7^ 0) 

for the variable elimination 



procedure for a non-interacting set Sa- Here "only" indicates up to multi- 
plication by a constant. 
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of denominators. The denominators can be chosen to be S-positive and multiphcation 
by the denominators does not change the positivity of sohitions. Fix a maximal set of 
n independent combinations ^\ci, . . . ,Cs) = Yli=i^i'^i providing conservation laws that 
includes those corresponding to the full connected components of (that is, to cuts). 
For given total amounts, the steady state equations are complemented with the equations 
uj^ = ^'(ci, . . . , Cs), I = 1, . . . ,n. If the conservation law corresponds to a cut, then the elim- 
ination procedure ensures that CK^ii^a)^ • • • ; ^m{Ca)^ Cm+i, • • • , c^) = and the equation 
becomes redundant. 

Theorem 9.1. Consider a CRN with a non-interacting set Sa- Assume that Sa = S^U 
• • • U 5!I is a partition of Sa into disjoint sets such that G^j is connected and G^j admits 

a spanning tree for all j. If is a cut, assume that the spanning tree is rooted at some 
Si € Sa and otherwise assume that it is rooted at *. Let total amounts a;' be given for the 
n conservation laws. 

The non-negative steady states with positive values in Ca are in one-to-one correspon- 
dence with the positive solutions to 

<^>4C'a) =0, OJ' = e'(C^) := C'iMC'a), . . . , V^miC'a), C„+l, . . • , C,) 

for u = m -\- 1, . . . , s and I = I, . . . ,n. 

Proof. We have shown that any non-negative steady state solution with positive values 
for a G Ca must satisfy these equations. For the reverse, we apply the following to 
each connected component of Gs^. Consider a positive solution c = (cm+i, • • • ,Cs) to the 
equations ^u{Ca) = and = C'(C'a). For i = 1, . . . ,m, define Cj through Proposition 
or 



8.11 



8.9 



depending on whether Si belongs to a cut S^ or not. For positive rate constants 
and positive total amounts, q is non-negative (because the rational functions defining it 
are S-positive). By construction this procedure automatically ensures that conservation 



laws corresponding to cuts are fulfilled. Using Propositions 8.9 and 8.11 the values c, 



satisfy (4.2). Since <1>„(C^) = is the steady state equation Cu = after substitution of the 
eliminated variables, this equation is also satisfied and the same reasoning applies to the 
equation uj^ = ^'(ci, . . . , Cm)- Thus, (ci, . . . , Cm) is a solution to the steady state equations 
and satisfies the conservation laws corresponding to the total amounts loK □ 



This theorem together with Proposition 8.12 give the following corollary. 



Corollary 9.2. With the conditions of Theorem 9.1, assume further that each graph G^j 
is strongly connected. Then, the positive steady states of the system are in one-to-one 
correspondence with the positive solutions to <I>u(C^) = and uj^ = ^'(C^) for u = m -\- 
1, . . . ,s and I = 1, . . . ,n. Further, if a steady state solution satisfies Cj = for Ci £ Ca, 
then there exists some cj £ Ca such that cj = 0. 

In the main example, the set Sa = {Si, S4, S5, Sq, Sg, Sg} is the largest non-interacting 
subset of S and thus provides the maximal number of linearly eliminated concentrations. 
The initial steady state system of equations is reduced to three equations: For instance 
the one corresponding to cy = 0, and the two conservation laws oJ'^ = C2 + C4 + cq + ct -\- cg 
and uj = C3 + C5 + cq + C7 + Cg (which corresponds to — w^). Because of the conservation 
laws, the equations C2 = and C3 = are redundant. The elimination from cuts provides 
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Cjfc = o-fc/(Ti for k = 4,5,6, ci = a;Vi/(<Ti+<T4+<T5+<T6), cg = ag/as and cg = uPas/{as+ao) 
with (Tj S-positive polynomials in 02,03,07 and coefficients in ]R(ConU{aJ^,cj^}) for all i. 
The steady state equations are thus reduced to: 

= A;9(76((78 + (To)as - kiQuf'alaiC'j + A;ii(7icr9((j8 + ag) 

3 

u (Ticrg = a\a^C2 + 0"4Cr8 + (Jecrs + (JicrgCy + (Ticrg 
ujaiG^ = aiascs + (TsCTs + crecs + cncgCT + cricrg. 
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